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Abstract 


In this paper, an analytical method called the parametric iteration method 
(PIM) is presented for solving the second-order singular IVPs of Lane-Emden 
type, and its local convergence is discussed. Since it is often useful to have 
an approximate analytical solution to describe the Lane-Emden type equa- 
tions, especially for ones where the closed-form solutions do not exist at all, 
therefore, an effective improvement of the PIM is further proposed that is ca- 
pable of obtaining an approximate analytical solution. The improved PIM is 
finally treated as an algorithm in a sequence of intervals for finding accurate 
approximate solutions of the nonlinear Lane-Emden type equations. Also, 
we show how to identify an approximate optimal value of the convergence 
accelerating parameter within the frame of the method. Some examples are 
given to demonstrate the efficiency and accuracy of the proposed method. 


Keywords: Piecewise-truncated parametric iteration method; Truncated 


parametric iteration method; Parametric iteration method; Nonlinear 
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1 Introduction 


Lane- 


Recently, a lot of attention has been focused on the study of singular initial 
value problems (IVPs) in the second-order ordinary differential equations 


(ODEs). Many problems in mathematical physics and astrophysics c 


an be 


modelled by the so-called IVPs of the Lane-Emden type equation [2, 4, 14]: 
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where a and 0 are constants, f(x,y) is a continuous real valued function, and 
g(x) € C[0,co]. When f(x,y) = K(y) , g(x) = 0, Equation(1) reduces to 
the classical Lane-Emden equation which, with specified K(y), was used to 
model several phenomena in mathematical physics and astrophysics such as 
the theory of stellar structure, the thermal behavior of a spherical cloud of 
gas, isothermal gas sphere and theory of thermionic currents [2,4, 14]. 

Since, the Lane-Emden type equations have significant applications in 
many fields of scientific and technical world, a variety of forms of f(z, y) 
and g(x) have been investigated by many researchers (e.g., [3,16,17]). A 
discussion of the formulation of these models and the physical structure of 
the solutions can be found in the literature. Though the numerical solution of 
the Lane-Emden Equation (1), as well as other various linear and nonlinear 
singular IVPs in quantum mechanics and astrophysics [9], is numerically 
challenging because of the singularity behavior at the origin « = 0, but 
analytical solutions are much needed for physical understanding. Recently, 
many analytical methods were used to solve the Lane-Emden equation [7,8, 
10,18]. Those methods are based on either series solutions or perturbation 
techniques [1, 11-13]. However, the convergence region of the corresponding 
results is rather small. 

The strategy that will be pursued in this work rests mainly on establishing 
useful algorithms based on the parametric iteration method (PIM) [5,6, 15] 
for finding highly accurate solution of the Lane-Emden type equations that 
they 


e Overcome the main difficulty arising in the singularity of the equation 
at x = 0. 


e Provide us with a convenient way to modify the convergence region and 
rate of the solution. 


e Are simple to implement, accurate when applied to the Lane-Emden 
type equations and avoid tedious computational works. 


The examples analyzed in the present paper reveal that the newly developed 
algorithms are easy, effective and accurate to solve the singular [VPs of Lane- 
Emden type equation. 


2 Analysis of methods 


In this section, the PIM is described for solving Equation(1). This method 
provides the solution of Equation(1) as a sequence of iterations. The method 
gives rapidly convergent successive approximations of the exact solution if 
such a solution exists, otherwise approximations can be used for numerical 
purposes. The idea of the PIM is very simple and straightforward. To explain 
the basic idea of the PIM, we first consider Equation(1) as follows: 
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Lly(x)] + N[y(x)] = 9(2), (2) 
with 
Ly(2)] = y" (a) + =y'@), Nly(2)] = Fle), (3) 


where L denotes the so-called auxiliary linear operator with respect to y and 
N is a nonlinear operator with respect to y. The basic character of the PIM 
is to construct a family of iterative processes for Equation(1) as follows [15]: 
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where yo() is the initial guess and the subscript n denotes the n-th iteration, 
and h # 0 denotes the so-called auxiliary parameter which can be identified 
easily and efficiently by the technique proposed in this paper. Accordingly, 
the successive approximations y,(z) , n > 0 of the PIM in the auxiliary 
parameter will be readily obtained by selecting the initial approximation. 
Consequently, the exact solution can be obtained by using 

y(x) = lim yn(x). (5) 

noo 

It is interesting to note that for the linear Lane-Emden type equations, its 
exact solution can be obtained easily by only one iteration step due to the 
fact that the multiplier can be suitably identified, as will be shown in this 
paper later. 
Now we will have the following proposition for the iteration formula (4). 


Proposition 1. [f y(x) € C?[0,T], then, forx <T , 


[ (/ = | { (i) + 2y oka = y(n) ~ y(0). (6) 


Proof. Simple integration by parts. 


In the light of (4) and (6), we will have the following simple iteration 
formula: 


t2 


tne) = (1+ nC) — hao + f (/ - ") [f(t yn(t)) — g(@)]dt, (7) 


where yo = y(0). 
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Note that the expressions (4) and (7) demonstrate that the variational 
iteration method (VIM) [7] is a special case of the PIM when h = —1. The 
fact that the PIM solves the Lane-Emden type equations without correction 
functional and restricted variation can be considered as an advantage of this 
method over the VIM. 


The PIM (7) makes a recurrence sequence {y,(z)}. Obviously, the limit 
of this sequence will be the solution of Equation (1) if this sequence is con- 
vergent. 

In order to prove the sequence {y,(x)} is convergent, we construct a series 


yo(x) + [y1(@) — yo(@)| +--+ + [yn(%) — yn—1(@)| ++ (8) 


Noticing that 


Snti = yo(@) + [y(x) — yo(a)] +--+ + [yn(@) — Yn-1(#)] = yn(a), (9) 
the sequence {y,,(x)} will be convergent if the series is convergent. 
Theorem 1. Jf Niy] = f(x,y) is Lipschitz-continuous in [0,T] and 


g(x) € C[0,T] then the series of (8) is convergent, i.e., the sequence {Yn(x)} 
is convergent for x € [0,T). 


Proof. According to (7), note that 


h i : (+ ©) extn) ae <|h|MNe, (10) 


lyi(x) — yo(x)| = 


where 


M= 


= MAX a N IMbo)]— 9)! GY 


= max 
0<t<2<T 


From (7) and (10), and the assumption that |N[yn]—N[yn-1]| < Llyn—Yyn-1| 
where L denotes the Lipschitz constant of N[y(«)], it follows that 


ly2(x) — yi(@)| = [1 + Ally (x) — yo(%)| + |h| ML a ly1(t) — yo(t)|dt 
(\h| M La) ee 


u! 2! (2) 


Nee ti (\h| MLa)*+1 
ee 1+-Al-* 
Th) | (k+1)! ” 


IA 
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In view of (14), the convergence of the series (8) can be concluded for the 
solution domain x < T and |1+h| < 1 with the help of the mathematical 
software such as Maple. Therefore the series of (8) is absolute convergence, 
ie., the sequence {yp(x)} is convergent for x € [0,7]. 


3 A piecewise-truncated PIM 


The successive iterations of the PIM may be very complex, so that the re- 
sulting integrals in the relation (4) may not be performed analytically. Also, 
the implementation of the PIM generally leads to calculation of unneeded 
terms, which more time is consumed in repeated calculations for series solu- 
tions. Here, an effective modification of the PIM is introduced to eliminate 
these repeated calculations. To completely stop these repeats in each step, 
provided that the integrand of (4) in each of iterations is expanded in mul- 
tivariate Taylor series around x = 0. We propose the following improvement 
of the PIM (4), which is called the truncated PIM (TP): 


vaial#) = yale) +h F,(2,t)dt, n>0, (15) 
0 
where 
2 ww a 
(: : ‘) im () + =u) + Fe yn() - aco} = Fy (0,t) + O(a") +O"). (16) 


It is noteworthy to point out that the TP formula (15) can cancel all the 
repeated calculations and terms that are not needed as will be shown below. 
Furthermore, it can reduce the size of calculations. Most importantly, how- 
ever it is the fact that the TP algorithm (15) solves a Lane-Emden differential 
equation exactly if its solution is an algebraic polynomial up to some degree. 

In general, by using the TP formula (15), we obtain a series solution, 
which in practice is a truncated series solution. This series solution gives a 
good approximation to the exact solution in a small region of z. An easy 
and reliable way of ensuring validity of the approximations (15) for large z, 
ie. I; = (xj, 2j41] where Av = 2j4, — 2;, 1=0,1,2,---,N—1, with ap =0 
and ay = T. According to the relation (15), therefore, we can construct 
the following piecewise TP approximations (PTP) in the subintervals J;. On 
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[xo, 1], let 


Urvsni(e) = UiiL) + Riles F, m(z,t)dt, m= 0, 1, ao SD Ly rE [xo, x1], 


Fy m(2, t) + O[(@ — 20)™*1] + O[(t — 0)". 
(17) 


Then one can obtain the n,-order approximation yj n, (x) on [xo, 71]. 
On [21, x2] , let 


Y2,m41(2) = Y2,m(2) + hfe Fo m (a, t)dt, Mm = 0, 1, se lg — 1, Le [x1, Xa], 


y2,o(@) = Y1yny (21) + Yin, (@1)(@ — @1) = 1 + ey (@ — 21), 
(/ = ®) {al + 245 m(t) + f(t y2,m(t)) — 7) = Fo oh (Gy) 


+O[(x2 — 21)™*1] + O[(t — 21)™*1). 
(18) 


Also the ng-order approximation y2 n(x) on [%1, 2] can be obtained 
In a similar way, on [x;,2;41], 7 = 2,--- ,N —1, let 


Yittm+1 (2) = Yitajm(2) +h fe Fistm(a,t)dt, m=0,1,--- niga —1, 
LE [we tials 


Yi41,0(2) = Yin, (wi) + Visine (2-2) =q+ Co(x = Z;), 


(/ i :) {al + 2:4 1m(t) + Ft yiram(t)) — ao} = Fittm(2,t) 


Hence this implies the n;41-order approximation yj+41,n,,,(%) on [#;, 7541]. 
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Therefore, according to (17)-(19), the approximation of Equation (1) on 
the entire interval [0,7] will be calculated. It should be emphasized that 
the PIM and TF algorithms provide analytical solutions on [0,7], while the 
PTP technique provides analytical solutions in [x;,2;41], which are continu- 
ous at the end points of each interval, i-e., yin, (@i) = Ci = Yett nig, (i) and 


Yin, (Xi) = C; = Yitl niga (x3), t= 1, 2, ie ,N —1. 


Remark 1. In the case of failure of convergence of the PIM, the presence 
of the parameter could play a very important role in the frame of the PIM. 
Although, we can find a valid region for every physical problem by plotting 
the solution or its derivatives versus the parameter in some points [6,15], but 
an approximate optimal value of the convergence accelerating parameter h 
can be determined at the order of approximation by residual error: 


Re s(h) = [uote + Nl[yn(t)] - g(t) Pat, (20) 


dRes(h) 


ii 


One can easily minimize (20) by imposing the requirement 


4 Implementations 


To give a clear overview of the content of this study, the several Lane-Emden 
type equations will be studied. These equations will be tested by the above- 
mentioned algorithms, which will ultimately show the usefulness and accuracy 
of these methods. Moreover, the numerical results indicate that the approach 
is easy to implement. All the results here are calculated by using the sym- 
bolic calculus software Maple 17. 


Eexample 1. As a first example, let consider the following linear, the non- 


homogeneous Lane-Emden equation, i.e., Equation (1) with f(x,y) = y and 
gaysae 3 


” 2 7; 
ee ae ee (21) 


subject to conditions 


y(0)=1, y (0) =0. (22) 


The PIM has a very simple approach. Its concepts begin with dividing the 
left hand (21) into two parts, i.e., the auxiliary linear operator L and the 
nonlinear operator N as: 
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uw 


Ly(2)] =y"(2) + <y'(@) +y(z) and Nias. (23) 


This will allow us to construct a family of iterative processes for Eq.(21) as 
follows [14]: 


maste) =unte) +h f° [Erie 0] {ano + Funto-Peth at. 24 


Using simple integration by parts, similar to Proposition 1, implies that 


_ sin(z) 


y(0), (25) 


ax 


[ E sin(t — o) fy" (+ 2!) + uldat = v2) 
0 [& 


In the light of (24) and (25), the following PIM is: 


sin(x) 


eC CTC ene h a |< sint ») (Peat, (26) 


where yo = y(0) and yo(x) = y(0) + y'(0)(x). According to (26), therefore 
the following approximations with starting the initial guess yo(x) = 1 are: 


(2+ 2h)x + hsin(x) — ha(3 + 3x + 27)e~* 


n(o) = _ : (27) 
aye (2+ 4h + 2h?)x + (2h + h?) sin(x) — ha(6x + 2x? + 6+ 3ha + hx? + 3h)e7 
aN 2x 
(28) 
which the exact solution of Equation(21) yields for h = —-1, ice, 
«(34 3x + 27)e~* — sin(x 
y(a) = SS (29) 


22 


Example 2. As other example, we consider the nonlinear and non-homogeneous 
Lane-Emden equation, i.e., Equation (1) with f(z, y) = y? and g(x) = 6+2° 
[10]: 


w 2 - 
y +—y ty? =6+2°, (30) 
subject to conditions 


y(0)=0, y (0) =0. (31) 


o] 
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Now, our aim is to solve the equation (30) by means of the TP algorithm 
(15). According to (15), it is easily to obtain the following approximations 
of the TP with starting the initial approximation yo(x) = 0: 


yo(x) = —h2?, (33) 
Yn(x) = [1—(1 +h)? 1]e?, n>3. (34) 


If we choose h = —1, then the TP algorithm yields the exact solution 
y(a) = 2°. (35) 


It is interesting to note that the TP algorithm (15) can solve a Lane-Emden 
type equation exactly if its solution is an algebraic polynomial up to some 
degree. 


Example 3. As final example, we consider the nonlinear, homogeneous 
Lane-Emden-type equation, i.e., Equation(1) with f(x,y) = 4(2e¥ +e?) and 
g(x) = 0: 


2 y 
y + “v + 4(2e¥ + e2) =0, (36) 
subject to conditions 


y(0)=0, y (0) =0. (37) 


Here, the aim to solve the equation (36) by means of the above-proposed 
methods. Since the integration of the nonlinear term 4(2eY + e?) in Equa- 
tion (36) is not easily evaluated, thus the PIM requires a large amount of 
computational work to obtain few iterations of the solution (we can replace 
the nonlinear term with a series of finite components). However, one can 
used the modified PIM method, i.e., the TP algorithm (15). According to 
(15), it is easily to obtain the following approximations of (36) with starting 
the initial approximation yo(«) = 0: 


yi(x) = (38) 
yo(a) = 2ha®, (39) 
y3(a) = He + 2h), (40) 
ya(x) = (6h + 6h? + 2h3)x? + (2h? + h3)a4 (41) 


and so on. To investigate the influence of h on the convergence of the solution 
obtained via the truncated PIM, here we plot the curves of Yoo(0 ) and yi) 
, as shown in Figure 1. According to these curves, it is easy to discover the 
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valid region h , which corresponds to the horizontal line segments. Now, in 
the light of (20), an approximate optimal value of h can be determined by 
the following residual error: 


y20(t) 
2 


Res(h) = : {ugo(t) + Zupo(t) +420 + eM) VPA, (42) 


By imposing the requirement dRes(h)/dh = 0 and solving the resulting equa- 
tion, we can obtain the approximate optimal value for h = —1. Fig. 2 shows 
a comparison of approximation obtained using the 20th-order TP algorithm 
for h = —1 with the exact solution of Eq.(36), i-e., y(z) = —2In(1+ 2”). 


-3.990 24.0105 
-3.995 24.0054 
y"(0) 4 y)(0) 24 
-4.005 23.9954 
-4.010+4 : + + 23.9904 : 
-2 15 } 0.5 0 16-14-12 } 08 -06 -04 


Oe 04 0.6 0.8 1 


“15 


Figure 2: Approximate solution for Example 3 using the TP algorithm where the dotted- 
line: the 20th-order TP algorithm when h = —1 and symbol: the exact solution 


As observed, the TP algorithm (15) in solving Equation (36) gives good 
approximations to the exact solution in a small region of xz. In order to 
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enlarge the convergence region and rate of the series solution, here we im- 
plement the PTP (19) proposed in the section 3. According to (19), taking 
N = 1000 and njiy = 3,41 = 5,4 = 0,1,---N — 1, we can obtain the 
approximations of (36) on [0,100 . Figure 3 shows the absolute errors (the 
differences between the approximate values and the exact values) of the PTP 
solution for nj41 = 3, nig, = 4, Av = 0.1 and the approximate optimal 
value h = —1. From Figure 3, it is easily to found that the present approxi- 
mations are efficient for a larger interval. 


0.000207, 
0.0000207, 


0.000184 
0.000164 
0.000144 0.000015 
0.000124 

E(x) 0.000104 oe) 0.000010- 

0.000084 

0.000064 


0.000005 5} 
0.000047, 


0.000027) 


Figure 3: Shows the absolute errors (Ex (#) = |texact(x) — up(x)|, k = 4,5) of the 
PTP solution with Ax = 0 for Example 3 where left: the Abs. Err. of the 4th-order 
PTP solution for the approximate optimal value h = —1 and right: the Abs. Err. of the 
5th-order PTP solution for the approximate optimal value h = —1 


In closing our analysis, we point out that three concreted modeling equa- 
tions of second-order singular [VPs of the Lane-Emden type equation were 
investigated by using the algorithms proposed. The obtained results have 
shown noteworthy performance. 


5 Conclusions 


Application of the methods based on the PIM presented in this paper to the 
three Lane-Emden type equations indicates that for the linear Lane-Emden 
type equations, its exact solution, if such a solution exists, can be obtained 
easily by only one iteration step due to the fact that the multiplier can 
be suitably identified, that the TP algorithm can solve a nonlinear Lane- 
Emden differential equation exactly if its solution is an algebraic polynomial 
up to some degree, and that for nonlinear Lane-Emden type equations can 
be useful in general. The numerical results demonstrate that the PIM is a 
useful analytic tool for solving the Lane-Emden type equations. 
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